rm(list=ls())


# Check that required packages are installed:
want = c("foreign", "ggplot2", "MASS", "reshape2", "plyr", "dplyr", "ggrepel", 
         "lme4", "coda", "relaimpo", "relimp", "arm", "boot", "entropy", "here")
have = want %in% rownames(installed.packages())
if ( any(!have) ) { install.packages( want[!have] ) }
# load packages
junk <- lapply(want, library, character.only = TRUE)
rm(have,want,junk)

options(scipen=999)

#=============================================================================================================
# Load Image
load("allmodels.RData")
#=============================================================================================================


# Function to extract length of range
range.p <- function(x){
  range(x, na.rm = T)[2] - range(x, na.rm = T)[1]
}

disp.lr <- ddply(subset(data, lrg_ch != 3 & center == 0  & all_same_lrp == 0), 
                 .(syslab, id, lrg_ch), summarise,
                 lr.d = range.p(lrp),
                 tot = length(!is.na(lrp)),
                 sameside = mean(sameside_ch),
                 otherside = mean(otherside_ch))
disp.lr <- disp.lr[disp.lr$tot > 1, ]
disp.lr$group <- NA
disp.lr[disp.lr$sameside == 1, ]$group <- "In-group Parties"
disp.lr[disp.lr$otherside == 1, ]$group <- "Out-group Parties"

disp.ag <- ddply(disp.lr,.(syslab, group), summarise,
                 mean.d = mean(lr.d),
                 sd.d = sd(lr.d),
                 up.d = quantile(boot(lr.d, sample.mean, R = 1000)$t, p = 0.975),
                 lo.d = quantile(boot(lr.d, sample.mean, R = 1000)$t, p = 0.025))
disp.ag <- merge(disp.ag, ratio.lm.boot, by = "syslab")

cor.test(disp.ag$ri_l2_m[disp.ag$group == "In-group Parties"],
         disp.ag$mean.d[disp.ag$group == "In-group Parties"])
cor.test(disp.ag$ri_l2_m[disp.ag$group == "Out-group Parties"],
         disp.ag$mean.d[disp.ag$group == "Out-group Parties"])

quartz(type = 'pdf', file = 'generalization_range.pdf',
       width=11,height=6, dpi=300)
ggplot(disp.ag, aes(x = ri_l2_m, y = mean.d, weigths = 1/sqrt(sd.d*ri_l2_sd))) +
  geom_errorbar(aes(ymin = lo.d, ymax = up.d), col = "gray75") +
  geom_errorbarh(aes(xmin = ri_l2_lo, xmax = ri_l2_hi), col = "gray75") +
  geom_smooth(col = "gray40", method = "lm", se = F) +
  geom_point(colour = "black") + 
  geom_text_repel(aes(label = syslab), size = 3) +
  facet_wrap(~group, ncol = 2) +
  theme_bw() +
  xlab("Importance ratio (\U03C9) \n(Gray lines are 95% bootstrapped confidence intervals)") + 
  ylab("Standard Deviation Party Placement\n(Gray lines are 95% bootstrapped confidence intervals)")
dev.off()

rm(disp.lr, disp.ag, range.p)
